### DIVIDE GENES INTO GENE FAMILIES part 1
### need to prepare the protein file 
# makeblastdb -in Chiifu_v3.1.pep.fa -dbtype prot -title chiifu -parse_seqids -hash_index -out chiifu
# blastp -db chiifu -query Chiifu_v3.1.pep.fa -out result -outfmt 6

library(dplyr)
library(Biostrings)

blast_result<-read.delim('result',header=F,sep='\t',stringsAsFactor=F)
names(blast_result)<-c('qseqid','sseqid','pident','length','mismatch','gapopen','qstat','qend','sstart','send','evalue','bitscore')

protein_seqs<-readAAStringSet('Chiifu_v3.1.pep.fa')
gene_list<-sapply(strsplit(names(protein_seqs),split=' '),'[',1)
gene_list=gsub(pattern='>',replacement='',x=gene_list)
protein_list<-sapply(strsplit(names(protein_seqs),split=' '),'[',1)

info<-data.frame(gene_list=gene_list,protein_list=protein_list)
row.names(info)<-info$protein_list

blast_result$qgene<-as.character(info[blast_result$qseqid,'gene_list'])
blast_result$sgene<-as.character(info[blast_result$sseqid,'gene_list'])

blast_result$pairs<- ifelse(blast_result$qgene > blast_result$sgene, paste(blast_result$qgene,blast_result$sgene,sep='_'),paste(blast_result$sgene,blast_result$qgene,sep='_'))
blast_result_filtered<-blast_result%>%filter(evalue<0.00001)%>%filter(bitscore>50)%>%filter(qgene != sgene) %>% group_by(pairs) %>% arrange(desc(bitscore)) %>% dplyr::slice(1)
write.table(blast_result_filtered,file='result_nangua_filtered',row.names=F,sep='\t',quote=F)